…r optimisations

Changes (an empirical seed-count convergence study found the per-timestep
multistart to be heavily over-sampled - see 'seed_study.py' and the README):

* Seed screening ('screen then polish'): the objective function is first
  evaluated once at every seed (costs roughly the same as fully optimising
  one or two seeds), then only the 'seed_screen_top_n' best seeds plus
  'seed_screen_uniform_n' uniformly spread seeds are fully optimised with
  NLopt. Defaults (16+16 out of 400 seeds) reproduce the full 400-seed
  multistart minimum to within ~0.5% cost at ~10x less CPU.
  Set seed_screen_top_n=None to restore the original behaviour.

* Cache the ObjectiveFunction per process per timestep (it was previously
  constructed - including loading rotation/trench/velocity files - once
  per seed).

* Skip the isochron interpolation ('isopolate') and the loading of the
  ridge/isochron/isocob files when fracture zones are disabled (they are
  disabled in all current configs; this work was previously done at every
  timestep and then unused).

* Safety cap ('nlopt_max_eval_safety', default 1000) on COBYLA objective
  evaluations: COBYLA occasionally fails to converge and can cycle for
  thousands of evaluations, stalling an entire MPI rank and hence the
  whole timestep.

* New resumable study/profiling harness 'seed_study.py' (prep/probe/run/
  report) to re-run the seed-count convergence study on any model.

* Environment variable overrides (OPTAPM_START_AGE, OPTAPM_END_AGE,
  OPTAPM_MODELS, OPTAPM_SCREEN_TOP_N, OPTAPM_SCREEN_UNIFORM_N,
  OPTAPM_SERIAL) for quick test runs without editing the config.

The optimisation is dominated by the net rotation (NR) minimisation
(subduction zone migration and continent speeds are not independent of it -
see Muller et al. 2022, Solid Earth), so a defensible uncertainty envelope
for the optimised reference frame consists of end-members in the NR bounds
rather than perturbations of the component weights:

  1. The no-net-rotation reference frame (zero-NR end-member) - already
     produced by every run.
  2. The 'best' run (default) - NR bounds (0.08, 0.20) deg/Myr as before.
  3. The 'nr_max' run - NR upper bound relaxed to 0.30 deg/Myr, the top of
     the range permitted by asthenospheric shear / seismic anisotropy
     constraints (Conrad & Behn 2010; Becker 2006).

The variant is selected via 'run_variant' in Optimised_config.py or the
OPTAPM_VARIANT environment variable, and (for variants other than 'best')
is appended to the model name so each variant writes its own output files.

…fault)

The original trench migration component drives trench-normal migration
velocities toward zero. Since zero net rotation also implies near-zero
trench migration, that constraint is largely redundant with the net
rotation minimisation. It is also at odds with observations: most trenches
roll back toward the ocean basin behind them rather than advancing toward
the overriding plate - across reference frames, 62-78% of trench segments
retreat, with mean trench-normal velocities of +1.3-1.5 cm/yr and medians
of +0.9-1.3 cm/yr (Schellart et al. 2008, Earth-Science Reviews;
Williams et al. 2015, EPSL).

The new default scheme ('trench_migration_scheme = "rollback"') drives
per-trench orthogonal velocities toward a target of +10 mm/yr retreat and
tightens the bounds on the mean orthogonal velocity from (-30, 30) to
(0, 20) mm/yr - the global mean must be retreating, but at less than
2 cm/yr. Individual trenches can still advance; only the mean is required
to retreat. This makes trench kinematics an independent constraint that
competes with the net rotation minimisation (instead of duplicating it).

Set 'trench_migration_scheme = "minimise"' (or OPTAPM_TM_SCHEME=minimise)
to restore the original behaviour.

…ion weight

Relaxing the net rotation bounds alone is not sufficient to produce a
directional end-member: the bounds are hard penalty walls that only affect
timesteps where the optimum presses against the 0.20 deg/Myr ceiling
(mostly in deep time), whereas inside the walls the solution is set by the
smooth trade-off between the NR cost and the other components. Halving the
NR weight (inverse weight 2.0) shifts that trade-off at every timestep,
letting the trench rollback and plate velocity constraints pull harder
against the net rotation minimisation.

Single-timestep tests (5 Ma) confirm: bounds-only relaxation merely
selected a different near-degenerate minimum (NR 0.123 deg/Myr), whereas
bounds + halved weight produced the intended systematically
stronger-rollback, higher-NR end-member (NR 0.177 deg/Myr, trench retreat
strengthening from +1.8 to +3.2 mm/yr, 68% of segments retreating).

…ough time)

At the end of every run (config flag 'generate_diagnostics', default True)
the workflow now computes per-timestep summary statistics of the optimised
model and writes them to 'model_output/':

* <model_name>_diagnostics.csv - per timestep: median and MAD of the net
  rotation rate (deg/Myr, sampled at 1 Myr sub-steps within each timestep),
  and mean/median/MAD of the trench-normal migration velocities (mm/yr,
  positive = retreat) across all subduction zone segments, plus the
  percentage of segments retreating.
* <model_name>_net_rotation.png - NR (median +/- MAD) through time, with
  the 0.20 deg/Myr preferred upper bound and the 0.26 deg/Myr
  Conrad & Behn (2010) limit marked.
* <model_name>_trench_migration.png - trench-normal migration
  (mean +/- MAD) through time, with the observed present-day mean retreat
  range (+0.9 to +1.5 cm/yr, Schellart et al. 2008) marked.

The same statistics/plots are also generated for the no-net-rotation model
(outputs '<model_name>_NNR_*') - showing subduction zone kinematics in an
NNR world (the zero-NR end-member of the uncertainty envelope), with the
NNR net rotation plot doubling as a sanity check (~zero at all times).

Diagnostics can be (re)generated standalone on an existing optimised model
via 'python model_diagnostics.py'.

Also added a 'References for parameter choices' section to the README and
expanded the justifications (with references and validation tables) for
all recently introduced parameters: seed screening counts, COBYLA
evaluation cap, net rotation bounds and weights per run variant, and the
trench rollback target and bounds.